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We obtain in TM polarization an analytical expression of the scattering ma- 
trix of one infinitely conducting metallic lamellar grating with subwavelength 
slits. The theory is based on the Monomode Modal Method which consists 
in considering only one propagative mode in grating slits. Two expressions 
are exposed. The first one comes directly from the theory equations and the 
second one clearly reveals the Airy-like form of the scattering matrix terms. 
The theory is validated on a multi-grating object and the stability of the 
numerical results are shown at the same time. This work provides a basic 
and very efficient theoretical tool to calculate the diffraction by a stack of 
subwavelength metallic gratings. ©2013 Optical Society of America 
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1. Introduction 

Since the discovery of the extraordinary optical transmission through subwavelength metal- 
lic gratings (SMG) [1,2], a lot of numerical methods [3-6] have been applied to understand 
this singular physical phenomenon and to calculate the transmitted intensities. Now, it is 
known that such periodic sets of plasmonic or cavity resonators behave as Fabry-Perot in- 
terferometers. Today, the scientist's interest turns widely to applications, such as tunable 
optical filters [7-9] or micro-polarizer exhibiting optical activity [10,11]. Some of these opti- 
cal properties require complex structures made of a stack of SMG [12-15] which often implies 
tricky and time-consuming computations. Thus, from fundamental physical considerations, 
it seems important to enhance the theoretical analysis in order to obtain very efficient com- 
puting codes to increase understanding of such structures and to easily discuss the influence 
of opto-geometrical parameters. 
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For planar periodic objects like gratings surrounded by two semi-infinite homogeneous 
regions, the common and well-known theories such as the Classical Modal Method (CMM) 
[16,17], the Rigorous Coupled- Wave Method (RCWM) [3,18,19] or the Differential method 
(DM) [20,21], usually lead to the calculation of two kinds of matrices relying the fields 
in homogeneous regions: the Scattering matrix (S-matrix) or the Transmission matrix (T- 
matrix). From them, iterative processes like the Scattering Matrix Propagation Algorithm 
(S-algorithm) [21,22] are then used to compute the diffracted fields in multi-layer devices. For 
the case of gratings made of infinitely conducting metal, the CMM has remained for the last 
decades the most relevant method to solve the diffracting problem. It consists in considering 
the grating slits as waveguides and to express the cavity fields as a combination of waveguide 
modes. The initial formalism of the theory imposes to take as many modes (propagative or 
evanescent) in grating cavities as Rayleigh terms in homogeneous regions. But the transverse 
fields of evanescent modes make the method unstable (one matrix to invert becomes singular) 
when the slit width decreases and the number of harmonic terms increases. These numerical 
problems have limited its applications for a long time. By relevant matrix transformation, 
Gralak et al have obtained non singular system to invert [23] . The most studied configuration 
is when the cavity width is small enough to have only one propagative mode. The incident 
wavelength being larger than the cut-off of all the modes, except the fundamental one. 
When only one mode propagates in slits as in studied monoperiodic gratings lighted in TM 
polarization (TMq fundamental mode), this method has been named the Monomode Modal 
Method (MMM). This theory has been used before so far only to analytically express the 
propagative transmitted or reflected diffracted fields from which the diffracted efficiencies, 
transmittance and reflectance are deduced [24-26]. 

We propose to extend the MMM to the analytical calculation of the SMG S-matrix. It is 
worth noticing that the T-matrix cannot be analytically obtained since numerical inversion 
of one large matrix is required [12]. We only consider monoperiodic gratings lighted by a 
TM planewave. Therefore, the theory is restricted to thick gratings and wavelengths larger 
than the cut-off of the second cavity mode (TMi). At first, we give the analytical expression 
of the SMG S-matrix directly deduced from the MMM equations. After that, we propose a 
more physical writing of the S-matrix terms from an Airy-like formulation which includes 
Fresnel coefficients at grating interfaces. Finally, the method is validated on a multi-grating 
device with the help of the S-algorithm and we verify that the numerical results are without 
numerical divergences. 

2. Scattering matrix of one subwavelength metallic grating 

We establish in this section the analytical expression of the S-matrix for one SMG sur- 
rounded by two semi-infinite homogeneous regions (see fig. 1). We assume that the metal is 
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infinitly conducting and that refractive indices are noted n 2 for grating cavities and rij for 
homogeneous regions j G {1,3}. The grating period is denoted d, its thickness h = — z\ 
and the cavity width w. In TM polarization, we know that the transverse fields are written 
as Fourier- Ray leigh (FR) expansions in homogeneous regions (j): 

-iuit 



= ej—— e iapX Afj^^-^ + £C7') e -*# J (*-*i) , j e {1, 3} (I) 



with ^ = [ejfi u 2 — «p] and a p = a + p2n/d. 

We note that the wavelength validity domain is [A c , +oo], where A c = 2w is the cut-off 
of the TM\ mode and that the TMq mode has no cut-off. In fact, the cavity resonances 
related on high order modes (TM m , m G [2,+oo]) cannot be excited since their resonance 
wavelengths are widely distant from this wavelength range. In the same way, their cut-off 
wavelengths are also widely smaller than the considered wavelengths since the grating is 
made of subwavelength cavities (w « A). Consequently, we can suppose that propagative 
and contra-propagative modes totally vanish during their evanescent propagation (specially 
for high /i-values) and cannot establish stationnary resonances. For all these reasons, it seem 
judicious to consider in our theoretical analysis only the unique propagative mode TM . 

Consequently, the transverse field in cavities is written as the sum of propagative and 
contra-propagative TM Q modes with amplitudes A and B: 
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with ko = 2tt/\. We emphasize that the choice of the phase origin in field expressions (1) and 
(2) is crucial to obtain an analytical expression of the S-matrix only depending on h. The 
projection of the field continuity relations for both electric and magnetic fields at grating 
interfaces on FR-orders and cavity mode transverse fields leads to 4 equation sets [25]: 

Af + Bf> = (lu + B) g p , \/p G Z (3) 
4° + B p ] = { I+ ® u ) 9p, VP^ Z ( 4 ) 
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E ^? [4 3) - B f ] ] 9l = n* (An - B) (6) 
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where u = expiik^n^K) with ko = 2tt/X, and the exponent * denotes the complex conjugate. 
The overlap integrals g p between the cavity modes and FR-orders are: 

9p = V7 smc ( a ff) ( 7 ) 
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The quantity rjp = koCj/^p is the admittance of FR p-order. 

Injecting Eq. (3) into (6), then (4) into (5) leads to obtaining an equation set linking the 
mode amplitudes and the incident FR-orders. It writes in matricial form as: 
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cm . 



n 2 
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where C^' are very important effective coefficients related to a sum of overlap integrals 
between the TM mode and FR modes ponderated by the FR admittance: 



(9) 



As usual, the FR-expansions are in the following truncated into 2iV + 1 harmonic terms 
in order to obtain matrices with finite size. Introducing column matrices [A^')] and [J3^] 
containing FR-amplitudes A p ^ and Bp^ 1 respectively, the analytical inversion of the system 
(8) which only requires inversion of a 2 x 2 matrix, leads to simple relations relating A and 
B to all the Ap 1 ^ and B^ amplitudes: 
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where V is a matrix of size 2 x (2N + 1). The elements of its four blocks can be explicitely 
written: 
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with p G [— N, N]. Finally, an analytical expression of the S-matrix is basically obtained 
from Eqs. (3) and (4) by simply replacing A and B by their expressions given by Eq. (10) 
and (11) to (14): 

/ r 4(3)1 \ / r 4(di \ 
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with (p, q) G [-N, N] x [-N, N] and S pq the Kronecker symbol. The size of S-matrix is 
2(27V+ 1) x 2(2JV + 1). We notice that (^ 22 ) pa = (Sn) p Jv^- 



At this stage of our calculation, a rapid verification can be obtained. If only one incident 



FR-order is sent on the structure: A 
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0, in this case the diffracted 



amplitudes are the reflection and transmission coefficients: A p = t p and B p 



r p . Their 
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analytical expression are given by the relation t p = (Su) and r p = (S21) p0 which coincide 
with published ones [25]. In the same way, A = (Vu) and B = (P 2 i) provide the expresion 
of the cavity mode amplitudes. 



3. Another expression of S-matrix 

We can obtain a more physical expression of the S-matrix which lets clearly appear the 
reflection and transmission coefficients at interfaces. At first, we have to explicit the Fresnel 
coefficients at the interfaces between a homogeneous region (1) and a metallic waveguide 
(2) as shown in fig. 2. The fields are expressed as in eq. (1) and (2) but Z\ = z 3 = z', the 
interface coordinate on z-axis. We introduce the S-matrix of the interface: 




£(1.2) r (2,l) 
r (l>2) £(2,1) 




(20) 



Noticing N the truncation order of FR-expansions, we highlight that tM 2 ^ is a line 1 x (2A + 
l)-block, rM 2 ) is a square (2N + 1) x (2N + l)-block, rM 1 ) i s a scalar block and fM 1 ) i s a 
column (2 A + 1) x 1-block. We can easily prove that their elements are expressed by: 
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with p G [-N, N] and g G [-N, N]. We see that ip 1 ' 2 ^ and r'- 2 ' 1 '' are defined as usual amplitude 
ratios: t^' 2 ^ = and r^ 2,1 ) = i without incident g-order. But tp 2,1 ^ and r^ q 2 ^ are partial 

(1 2) 

reflection and transmission coefficients. For instance, r}, >q expresses the reflected g-order 
when only the incident p-order falls on the gratings. A summation is necessary otherwise 
(several incident p-orders). 

The S-matrix of one metallic grating is obtained by applying twice the S-algorithm be- 
tween the two S-matrices of each interfaces and the S-matrix related to the cavity mode 
propagation. We repeat that it is an iterative algorithm which consists in calculating de 
S-matrix of two adjacent layers according to their S-matrices S^ a ' and [22,27]: 
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where Id denotes the identity matrix. Before applying the S-algorithm, we remark that the 
phase origin of the cavity fields chosen in eq. (2) is not the same as the one used to explicit 
the S-matrix of one interface in eq. (20). Thus, we have to replace B by Bu in eq. (20) for 
the matrix at z — Z\. In the same way, A is changed to Au for the matrix at z = z 3 . 
The mode propagation S-matrix in cavities is obviously equal to uld- The S-algorithm 
is first applied on S'M and which leads to the matrix S^ 1 ' 2 ': 

( f (l,2) 2 (2,1) \ 

s ™ - { % ) « 3 °> 

A second application of S-algorithm on S'M 2 ) and S'M leads to the S-matrix S^ 1 ' 3 -* = S of the 
grating. Its blocks have a simple expression: 

4 2 ' 3) t (1 ' 2) « 

( S U) M = 1 _ r ( 2 ,l) r (2,3) M 2 ( 31 ) 
( Sl2 >P,q ~ I _ r (2,l) r (2,3) M 2 °P« ( 32 ) 
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P>1 \ _ r (2,l) r (2,3) u 2 

This formulation is equivalent to the one given in eqs. (16) to (19) but clearly reveals that 
the S-matrix terms can be expressed as Airy-like formulae. We also prove that the P-matrix 
can be written: 

t {1 ' 2) 

(V^P = l _ r (2 P l) r (2,3) n 2 (35) 

y Vl2 >P ~ l _ r (2,l) r (2,3) M 2 ( 36 ) 
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C^L = Z j,, U , :!L o (38) 
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4. A convergence test 

Before validating the theory on multilayer devices, we have analyzed the numerical stability 
of the analytical S-matrix computation on a single SMG. Many numerical methods based 
on similar formalisms to the CMM, the RCWM or the DM are fundamentally unstable. 
Precisely, they often diverge when the truncation order of FR-expansions N and/or the 
grating thickness increase. The S-algorithm and the Fast Fourier Factorization [20,21] bring 
solutions to these numerical problems. In order to evaluate the stability of our method, we 
have to define a convergence criteria. For every value of N and h, the computed quantity 
R + T is exactly equal to 1 (energy balance criteria) and so cannot be chosen to estimate 
convergence. Thus, we introduce the relative accuracy <r(N) = [T(N) — T(N max )]/T(N max ), 
between the transmittance value at N and the one at N max = 100 (N < N max ). The 
numerical analysis is done on the same structures studied in [12]. The figure 3 shows cr(N) 
according to N for different values of h and when the transmittance is evaluated to the first 
peak maxima (at d/X = 0.385 for h = 8/7). The transmittance versus d/X is also plotted 
in fig. 3 for every h value (plotted only on the wavelength range of the first peak for h = 2 
to 4). The other parameters are fixed to d — 1, w — 1/7, n\ = n% = 713 = 1 and Q{ nc = 0°. 
We see that the numerical results converge and remain stable whatever the h value. In 
fact, it is shown that the numerical instabilities in the classical modal method result from 
vanishing overlap integrals between high-order cavity modes and high-order FR-amplitudes, 
which leads to non-invertible matrices. Gralak et al [23] have proposed a solution to avoid 
such numerical problems by considerations on the matrix form. Stable formalism can be 
intrinsically obtained by considering only the propagating cavity mode as in our equations. 
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The convergence to zero of the overlap integral given by eq. (7) does not affect the stability of 
the analytical S-matrix terms: The coupling coefficients C)v) defined in eq. (9) and expressed 
from overlap integrals appear in the denominator of eqs. (16) to (19) and never tends to zero 
(but its modulus can tends to infinity at the Rayleigh wavelengths [25]). More generally, it is 
proven that the denominator in S-matrix terms vanishes only at resonance conditions when 
the free oscillation problem is solved (searching of complex resonance frequencies) [25]. 

5. A multilayered structure 

The numerical analysis of the diffraction by a multilayer device necessarily requires the 
calculation of its S-matrix. The use of the iterative S-algorithm (see eq.(25) to (29)) allows 
the computation of this global S-matrix of the studied structure from the S-matrices of 
each layer. The number of layers is denoted L. One layer can be a metallic grating (see 
fig. 1) or a homogeneous region as for the example of stacks of metallic gratings separated 
by homogeneous regions. Thus, it is necessary to explicit the S-matrix of a homogeneous 
layer. The fields in the homogeneous cavity are written as FR-expansions. The P-matrix of 
a homogeneous layer is a 4 x 4 diagonal block matrix. Its elements are expressed as in eqs. 
(35) to (38) but with u = u p and r^ 2 '^ = r p 2 , j G {1, 3}. All the reflection and transmission 
terms are the classical Fresnel coefficients for FR-orders: 

(1) _ (2) 

(1,2) _ VP VP fQH) 

'P (1) . (2) \° V > 

T)p +Vp 

l P (1) , (2) 

where (1) and (2) refer to the two homogeneous regions on both sides of one interface. The 
S'-matrix is also a 4 x 4 diagonal block matrix. Its elements are expressed as in eqs. (31) to 
(34) but with g p = 1, (Pk,i)p,q = (Pk,i) P 5 P ,q, (k,l) G {1,2} and with notations introduced in 
eqs. (39) and (40). 

We validate our theory on a stack of identical metallic gratings lighted in normal incidence. 
We consider the same grating as the previous one. The gratings are separated by identical 
homogeneous layers of thickness hhom = 4/7 and filled with air. We have studied four struc- 
tures from L = 1 to L = 7 where L = 2L g — 1 is the total number of layers and L g G N 
is the number of metallic gratings (there are L g — 1 homogeneous layers). N is fixed to 20. 
The transmittances are plotted in fig. 4. We obviously find again the transmittance curve 
for one and four gratings as studied in [12]. We also verify that new resonance peaks appear 
when a grating is added. In fact, some peaks are due to degeneracy splitting by coupling of 
two surrounded gratings, each of them behaving as a resonator. Other peaks are explained 
by the resonances of the homogeneous cavity created between two gratings and behaving as 
a Fabry-Perot resonator. 
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6. Conclusion 

We have obtained in TM polarization an analytical form for the S-matrix of a lamellar 
metallic grating because only one cavity mode is taken into account. In fact, the 2 x 2-matrix 
which links propagative and contra-propagative TM cavity modes can easily be inverted. 
As example, considering two propagative modes induces a tricky analytical inversion of a 
4 x 4-matrix. A generalization to a finite number of propagative slit modes leads to a semi- 
analytical formulation that we will present in a future work. We have also verified that the 
theory is numerically stable: the computations do not diverge when the number of Fourier- 
Rayleigh harmonics and the grating thickness increase. The main interest of the method 
remains the study of multilayer systems. It has been validated on a stack of metallic gratings. 

The theory may also be extended to the TE case but its wavelength validity domain is 
restricted to the wavelength interval for which one propagative mode exists, i.e. between TEq 
and TEi mode cut-off wavelengths. For the TM case, we restipulate that the TM mode 
has no cut-off which lends huge validity domain. 

In our future work, we shall apply this formalism to bi-periodic devices in order to develop 
nanostructures with polarization effects like, for example, an optical activity. We shall in 
addition present the semi-analytical calculation of the stable scattering matrix valid for 
metallic gratings with arbitrary cavity width including several propagative modes. 

Acknowledgements: I would like to thank Donna L'Hote from the Centre de linguistique 
appliquee (CLA) of Besancon (France) for her helpful advice. 
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List of Figure Captions 

Fig. 1. (Color online) Studied multilayer device: stack of metallic gratings. 

Fig. 2. (Color online) Schematic represantation of interface between an homogeneous region 
and a metallic waveguide for Fresnel coefficient calculation. Only the TM mode is considered 
in the waveguide. 

Fig. 3. (Color online) Convergence test of cr(T) according to the truncation Fourier- Rayleigh 
order N and for different values of h. cx(T) is evaluated at the first order resonance peak in 
shown transmittance curves versus A. 

Fig. 4. (Color online) Transmittance versus wavelength for multilayer devices. L denotes the 
number of layers. L = 1 for one grating ; L = 3 for two gratings and one homogeneous layer, 
etc. 
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Fig. 1. (Color online) Studied multilayer device: stack of metallic gratings. 
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Fig. 2. (Color online) Schematic represantation of interface between an ho- 
mogeneous region and a metallic waveguide for Fresnel coefficient calculation. 
Only the TM mode is considered in the waveguide. 
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Fig. 3. (Color online) Convergence test of er(T) according to the truncation 
Fourier- Ray leigh order N and for different values of h. cr(T) is evaluated at 
the first order resonance peak in shown transmittance curves versus A. 
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Fig. 4. (Color online) Transmittance versus wavelength for multilayer devices. 
L denotes the number of layers. L = 1 for one grating ; L = 3 for two gratings 
and one homogeneous layer, etc. 
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